home *** CD-ROM | disk | FTP | other *** search
/ HPAVC / HPAVC CD-ROM.iso / FRCOMP.ZIP / fracpack.c next >
C/C++ Source or Header  |  1992-11-12  |  13KB  |  327 lines
  1. Here is some C source for fractal image compression.
  2. Notes:
  3. 1)  You will need DJGPP (a C compiler for DOS ) to handle the arrays.
  4. 2)  It is available for anonymous ftp at
  5. ftp.math.niu.edu:/pub/msdos/djgpp.stuff/v109   (i think, look around)
  6. 3)  A 386+387 or higher is REALLY required.
  7. 4)  This should compile on other systems... graphics stuff will need to 
  8. change though...
  9. 5)  Images are 256x256x256 shades of gray.  Totals 65536 bytes. 1 byte/pixel
  10. 6)  I'm not a C programmer, so don't complain about the code or lack of
  11. documentation.  
  12. 7)  Get the paper /pub/inls-ucsd/fractal-2.0.tar from lyapunov.ucsd.edu
  13. It is a very good description of the process by yuval fisher.
  14. 8)  Decompression code will be sent in another message.
  15. 9)  Enjoy... there is also some later work of mine in /pub/young-fractal
  16. at the same anonymous ftp site.
  17. 10) It will be updated soon... (2+ weeks)
  18.  
  19.  
  20. ----------------- Cut Here -------------------
  21.  
  22. /*****************************************************************************
  23.                          FRACPACK - a program for
  24.                          FRACTAL IMAGE COMPRESSION
  25.                           image.krd ==> image.ifs
  26. *****************************************************************************/
  27.  
  28. #include <stdio.h>
  29. #include <stdlib.h>
  30. #include <graphics.h>
  31. #include <math.h>
  32. int main(int argc, char **argv)
  33. {  
  34.     unsigned char IMAGE [256][256], Range[8][8][8], Domain[8][8];
  35.     int REDGE[32][32], DEDGE[241][241];
  36.     FILE *in, *out;
  37.     char *inf, *outf, rmsstr[13];
  38.     int xsize=256, ysize=256, x, y, dx, dy, rx, ry, bestdx, bestdy, inrx, inry;
  39.     int class1, class2, class3, class4, i, besti, nflips=8, rsize, rsize2,
  40. plotdx, plotdy, usegraph;
  41.     long int suma,  sumb,  sumab,  suma2,  sumb2, sumx, sumx2;
  42.     long int DSUMS[241][241][2];
  43.     float   fsuma, fsumb, fsumab, fsuma2, fsumb2, fmagica;
  44.     float bests, besto, scale, offset, rms, bestrms, ifs_tolerance;
  45.     struct trans_t {
  46.       unsigned char dx;
  47.       unsigned char dy;
  48.       signed char scale;
  49.       short int offset : 13;
  50.       unsigned short int flip : 3;
  51.       } trans[1];
  52.  
  53.    
  54.     if (argc < 4)
  55.         {
  56.         printf("usage: fracpack rms infile.ext outfile.ext [graphics flag]");
  57.         return 1;
  58.         }
  59.     ifs_tolerance = atof(argv[1]);
  60.     inf = argv[2]; outf = argv[3];
  61.     if (argc == 5) usegraph = -1;
  62.     else usegraph = 0;
  63.     if ((in = fopen(inf, "rb")) == NULL)
  64.         {
  65.         fprintf(stderr, "Cannot open input file.\n");
  66.         return 1;
  67.         }
  68.     if ((out = fopen(outf, "wb")) == NULL)
  69.         {
  70.         fprintf(stderr, "Cannot open output file.\n");
  71.         return 1;
  72.         }
  73.     fclose(out);
  74.     printf("Compressing %s to %s with an rms of
  75. %4.2f\n",inf,outf,ifs_tolerance);
  76.  
  77.     ifs_tolerance=ifs_tolerance*ifs_tolerance;
  78.     if (usegraph)
  79.         {
  80.         GrSetMode(GR_default_graphics);
  81.         for (x = 0; x < xsize; x++)
  82.         GrSetColor(x,x,x,x);
  83.         }
  84.  
  85. /* Get .KRD bitmap 8-bit-grey-scale pixels, put in IMAGE array.  */
  86.     for (y = 0; y < ysize; y++)
  87.     for (x = 0; x < xsize; x++)
  88.         {
  89.         IMAGE[y][x] = fgetc(in);
  90.         if (usegraph)
  91.             {
  92.             GrPlot(x,y,IMAGE[y][x]);
  93.             GrPlot(x+300,y,IMAGE[y][x]);
  94.             }
  95.         }
  96.    fclose(in);
  97.  
  98. /*************************************************************************/
  99. /*                            Classify Range's                           */
  100. /*************************************************************************/
  101.     if (!usegraph) printf("Classifying RangeS\n");
  102.     for (ry = 0; ry < ysize; ry+=8)
  103.     for (rx = 0; rx < xsize; rx+=8)
  104.         {
  105.         for (y = 0; y < 8; y++)
  106.         for (x = 0; x < 8; x++)
  107.             Range[0][y][x] = IMAGE [ ry  +y ] [ rx  +x ];
  108.  
  109.         class1 = (Range[0][0][0]==Range[0][0][1]) +
  110. (Range[0][0][1]==Range[0][0][2])
  111.                + (Range[0][0][2]==Range[0][0][3]) +
  112. (Range[0][0][3]==Range[0][0][4])
  113.                + (Range[0][0][4]==Range[0][0][5]) +
  114. (Range[0][0][5]==Range[0][0][6])
  115.                + (Range[0][0][6]==Range[0][0][7]);
  116.  
  117.         class2 = (Range[0][7][0]==Range[0][7][1]) +
  118. (Range[0][7][1]==Range[0][7][2])
  119.                + (Range[0][7][2]==Range[0][7][3]) +
  120. (Range[0][7][3]==Range[0][7][4])
  121.                + (Range[0][7][4]==Range[0][7][5]) +
  122. (Range[0][7][5]==Range[0][7][6])
  123.                + (Range[0][7][6]==Range[0][7][7]);
  124.  
  125.         class3 = (Range[0][0][0]==Range[0][1][0]) +
  126. (Range[0][1][0]==Range[0][2][0])
  127.                + (Range[0][2][0]==Range[0][3][0]) +
  128. (Range[0][3][0]==Range[0][4][0])
  129.                + (Range[0][4][0]==Range[0][5][0]) +
  130. (Range[0][5][0]==Range[0][6][0])
  131.                + (Range[0][6][0]==Range[0][7][0]);
  132.  
  133.         class4 = (Range[0][0][7]==Range[0][1][7]) +
  134. (Range[0][1][7]==Range[0][2][7])
  135.                + (Range[0][2][7]==Range[0][3][7]) +
  136. (Range[0][3][7]==Range[0][4][7])
  137.                + (Range[0][4][7]==Range[0][5][7]) +
  138. (Range[0][5][7]==Range[0][6][7])
  139.                + (Range[0][6][7]==Range[0][7][7]);
  140.         class1 = class1*1000 + class2*100 + class3*10 + class4;
  141.         REDGE[ry>>3][rx>>3] = class1;
  142.         }
  143.  
  144.  
  145. /*************************************************************************/
  146. /*               Classify Domain's and compute DSUMS array               */
  147. /*************************************************************************/
  148.         if (!usegraph) printf("Classifying DomainS\n");
  149.         for (dy = 0; dy < ysize - 15; dy++)
  150.         for (dx = 0; dx < xsize - 15; dx++)
  151.             {
  152.             DSUMS[dy][dx][0] = DSUMS[dy][dx][1] = 0;
  153.             if (usegraph) GrPlot(dx+308,dy+8,384);
  154.             for (y = 0; y < 16; y+=2)
  155.             for (x = 0; x < 16; x+=2)
  156.                 {
  157.                 Domain[y>>1][x>>1] = (IMAGE[dy+y  ][dx+x  ]
  158.                                     + IMAGE[dy+y  ][dx+x+1]
  159.                                     + IMAGE[dy+y+1][dx+x  ]
  160.                                     + IMAGE[dy+y+1][dx+x+1]) >> 2;
  161.  
  162.                 DSUMS[dy][dx][0] += Domain[y>>1][x>>1];
  163.                 DSUMS[dy][dx][1] += Domain[y>>1][x>>1]*Domain[y>>1][x>>1];
  164.                 }
  165.  
  166.         class1 = (Domain[0][0]==Domain[0][1]) + (Domain[0][1]==Domain[0][2])
  167.                + (Domain[0][2]==Domain[0][3]) + (Domain[0][3]==Domain[0][4])
  168.                + (Domain[0][4]==Domain[0][5]) + (Domain[0][5]==Domain[0][6])
  169.                + (Domain[0][6]==Domain[0][7]);
  170.  
  171.         class2 = (Domain[7][0]==Domain[7][1]) + (Domain[7][1]==Domain[7][2])
  172.                + (Domain[7][2]==Domain[7][3]) + (Domain[7][3]==Domain[7][4])
  173.                + (Domain[7][4]==Domain[7][5]) + (Domain[7][5]==Domain[7][6])
  174.                + (Domain[7][6]==Domain[7][7]);
  175.  
  176.         class3 = (Domain[0][0]==Domain[1][0]) + (Domain[1][0]==Domain[2][0])
  177.                + (Domain[2][0]==Domain[3][0]) + (Domain[3][0]==Domain[4][0])
  178.                + (Domain[4][0]==Domain[5][0]) + (Domain[5][0]==Domain[6][0])
  179.                + (Domain[6][0]==Domain[7][0]);
  180.  
  181.         class4 = (Domain[0][7]==Domain[1][7]) + (Domain[1][7]==Domain[2][7])
  182.                + (Domain[2][7]==Domain[3][7]) + (Domain[3][7]==Domain[4][7])
  183.                + (Domain[4][7]==Domain[5][7]) + (Domain[5][7]==Domain[6][7])
  184.                + (Domain[6][7]==Domain[7][7]);
  185.  
  186.         class1 = class1*1000 + class2*100 + class3*10 + class4;
  187.         DEDGE[dy][dx] = class1;
  188.         if (usegraph) GrPlot(dx+308,dy+8,384);
  189.         }
  190.  
  191. /* Run through all non-overlapping 8x8 "R" blocks in the image  */
  192.     for (ry = 0; ry < ysize; ry+=8)
  193.     for (rx = 0; rx < xsize; rx+=8)
  194.         {
  195.         if (usegraph)
  196.             {
  197.             GrLine(rx  ,ry  ,rx  ,ry+7,384);
  198.             GrLine(rx  ,ry  ,rx+7,ry  ,384);
  199.             GrLine(rx+7,ry  ,rx+7,ry+7,384);
  200.             GrLine(rx  ,ry+7,rx+7,ry+7,384);
  201.             }
  202.         else printf("Entering Range Loop\n");
  203.         sumb = sumb2 = 0;
  204. /*************************************************************************/
  205. /*                   Grab Range   (and its "flips").                     */
  206. /*************************************************************************/
  207.         for (y = 0; y < 8; y++)
  208.         for (x = 0; x < 8; x++)
  209.             {
  210.             Range[0][y][x] = IMAGE [ ry  +y ] [ rx  +x ];
  211.             Range[1][y][x] = IMAGE [ ry+7-x ] [ rx  +y ];
  212.             Range[2][y][x] = IMAGE [ ry+7-y ] [ rx+7-x ];
  213.             Range[3][y][x] = IMAGE [ ry  +x ] [ rx+7-y ];
  214.  
  215.             Range[4][y][x] = IMAGE [ ry+7-y ] [ rx  +x ];
  216.             Range[5][y][x] = IMAGE [ ry+7-x ] [ rx+7-y ];
  217.             Range[6][y][x] = IMAGE [ ry  +y ] [ rx+7-x ];
  218.             Range[7][y][x] = IMAGE [ ry  +x ] [ rx  +y ];
  219.             sumb+=Range[0][y][x]; sumb2+=Range[0][y][x]*Range[0][y][x];
  220.             }
  221.  
  222. /*************************************************************************/
  223. /*                  Grab Domain, compare with Range(s)                   */
  224. /*************************************************************************/
  225.         bestrms = 10000000000.;
  226.         for (dy = 0; dy < ysize - 15; dy++)
  227.         for (dx = 0; dx < xsize - 15; dx++)
  228.             if ((DEDGE[dy][dx] == REDGE[ry>>3][rx>>3])||
  229.                 (((dx - 16 < rx)&&(rx < dx + 16))&&
  230.                 ((dy - 16 < ry)&&(ry < dy + 16))))
  231.                 {
  232.                 if (usegraph)
  233.                     {
  234.                     plotdy = dy; plotdx = dx;
  235.                     GrPlot(dx+308,dy+8,384);
  236.                     }
  237.                 for (y = 0; y < 16; y+=2)
  238.                 for (x = 0; x < 16; x+=2)
  239.                     Domain[y>>1][x>>1] = (IMAGE[dy+y  ][dx+x  ]
  240.                                         + IMAGE[dy+y  ][dx+x+1]
  241.                                         + IMAGE[dy+y+1][dx+x  ]
  242.                                         + IMAGE[dy+y+1][dx+x+1]) >> 2;
  243.  
  244.                 suma = DSUMS[dy][dx][0]; suma2 = DSUMS[dy][dx][1];
  245.                 fsuma=suma; fsuma2=suma2;
  246.                 fsumb=sumb; fsumb2=sumb2;
  247.                 fmagica = (float) (suma2 - suma*suma/64.);
  248.  
  249.                 for (i = 0; i < nflips; i++)
  250.                     {
  251.                     sumab = 0;
  252.                     for (y = 0; y < 8; y++)
  253.                     for (x = 0; x < 8; x++)
  254.                         sumab += Domain[y][x] * Range[i][y][x];
  255.                     fsumab = sumab;
  256.                     if (fmagica != 0.)
  257.                         scale = (fsumab - fsuma*fsumb/64.)/fmagica;
  258.                     else scale = 0;
  259.                     if (scale*scale < 1.44)
  260.                         {
  261.                         scale = (signed char) 127. * scale / 1.2;
  262.                         scale = 1.2 * scale / 127.;
  263.                         offset = (short int) (fsumb - scale*fsuma)/64.;
  264.                         rms = (fsumb2 + scale*(scale*fsuma2 - 2*fsumab +
  265. 2*offset*fsuma)
  266.                             + offset*(offset*64. - 2.*fsumb)) / 64.;
  267.                         if (rms > 10000) (i = nflips);
  268.                         if (rms < bestrms)
  269.                             {
  270.                             if (!usegraph) printf("%3i%4i%4i%4i%4i%2i
  271. %12.8f\n",REDGE[ry][rx],dx, dy, rx, ry, i, rms);
  272.                             besti = i;  bestrms = rms; bestdx = dx; bestdy =
  273. dy;
  274.                             bests = scale;  besto = offset;
  275.                             }
  276.                         if (rms < ifs_tolerance) goto gotbest;
  277.                         }
  278.                     }
  279.                 } /* End of Domain LOOP */
  280.  
  281. gotbest:
  282.         if (usegraph)
  283.             {
  284.             for (dy = 0; dy <= plotdy; dy++)
  285.             for (dx = 0; (dx <= plotdx) || ((dy < plotdy)&&(dx < xsize - 15));
  286. dx++)
  287.                 if ((DEDGE[dy][dx] == REDGE[ry>>3][rx>>3])||
  288.                 (((dx - 16 < rx)&&(rx < dx + 16))&&
  289.                 ((dy - 16 < ry)&&(ry < dy + 16))))
  290.                     GrPlot(dx+308,dy+8,384);
  291.             sprintf(rmsstr,"%8.4f",sqrt((double)bestrms));
  292.             GrTextXY(560,20,rmsstr,255,0);
  293.             }
  294.         if ((out = fopen(outf, "ab")) == NULL)
  295.             {
  296.             fprintf(stderr, "Cannot open output file.\n");
  297.             return 1;
  298.             }
  299.  
  300.         if (!usegraph)
  301.             {
  302.             printf("best rms=%12.8f,dx=%5i,dy=%5i,flip=%i for
  303.             rx=%5i,ry=%5i,s=%12.5f\n",  bestrms,
  304.             bestdx, bestdy, besti, rx, ry,bests);
  305.             }
  306.         trans[0].dx = bestdx;
  307.         trans[0].dy = bestdy;
  308.         trans[0].flip = besti;
  309.         trans[0].scale = 127. * bests / 1.2;
  310.         trans[0].offset = besto;
  311.         fwrite(trans, sizeof(struct trans_t), 1, out);  /* write out compressed
  312. file */
  313.         fclose(out);
  314.         if (usegraph)
  315.             {
  316.             GrLine(rx  ,ry  ,rx  ,ry+7,384);
  317.             GrLine(rx  ,ry  ,rx+7,ry  ,384);
  318.             GrLine(rx+7,ry  ,rx+7,ry+7,384);
  319.             GrLine(rx  ,ry+7,rx+7,ry+7,384);
  320.             }
  321.       }
  322.    if (usegraph) GrSetMode(GR_default_text);
  323. /* All done. Whew... */
  324.    return 0;
  325. }
  326.  
  327.